Please cite this using
Zhang, Z., McArdle, J. J., Wang, L., & Hamagami, F. (2008). A SAS interface for Bayesian analysis with WinBUGS. Structural Equation Modeling, 15(4), 705?C728.
TITLE 'Run WinBUGS from SAS: A Growth Model Example by Zhang & McArdle';
TITLE2 'Generate the Growth Model Data';
DATA Sim_LinGM;
* setting mathematical parameters;
MuL=10; MuS=5;
Sigma_e=1; Sigma_L=2; Sigma_S=1; rho=.5;
* setting statistical parameters;
N = 1000; seed = 20030930;
* need to setup arrays so we can have more variables;
ARRAY y_score{5} y1-y5;
* generating raw data;
DO _N_ = 1 TO N;
e_L=RANNOR(seed);
e_S=rho*e_L+SQRT(1-rho**2)*RANNOR(seed);
L_score=MuL+Sigma_L*e_L;
S_score=MuS+Sigma_S*e_S;
* now the indicator variables ;
DO t = 1 TO 5;
y_score{t} = L_score + t * S_score + Sigma_e*RANNOR(seed);
END;
KEEP y1-y5 L_score S_score;
OUTPUT;
END;
RUN;
PROC CORR DATA=Sim_LinGM;
RUN;
PROC MEANS DATA=Sim_LinGM;
VAR y1-y5 L_score S_score;
RUN;
DATA model;
INPUT model \$80.;
CARDS;/*Start of the model scripts*/
model
{#Model
for (i in 1:N){
LS[i,1:2]~dmnorm(Mu[1:2], Inv_cov[1:2,1:2])
for (t in 1:T){
y[i,t]~dnorm(MuY[i,t], Inv_sig2_e)
MuY[i,t]<-LS[i,1]+LS[i,2]*t
}
}
#Prior
Mu[1:2]~dmnorm(Mu0[1:2], Inv_cov0[1:2,1:2])
Mu0[1]<-0
Mu0[2]<-0
Inv_cov0[1,1]<-1.0E-6
Inv_cov0[2,2]<-1.0E-6
Inv_cov0[2,1]<-Inv_cov0[1,2]
Inv_cov0[1,2]<-0
Inv_cov[1:2,1:2]~dwish(R[1:2,1:2], 2)
R[1,1]<-1
R[2,2]<-1
R[2,1]<-R[1,2]
R[1,2]<-0
Inv_sig2_e~dgamma(.001,.001)
#Transform of the parameters
MuL<-Mu[1]
MuS<-Mu[2]
Cov[1:2,1:2]<-inverse(Inv_cov[1:2,1:2])
Sig2_L<-Cov[1,1]
Sig2_S<-Cov[2,2]
rho<-Cov[1,2]/sqrt(Cov[1,1]*Cov[2,2])
Sig2_e<-1/Inv_sig2_e
}
;
/*end of the model scripts*/
RUN;
DATA _NULL_;
SET model;
FILE 'C:\SASWinBUGS\GrowthModel.txt';
PUT model;
RUN;
/*Transform the SAS data file to the WinBUGS format
The SAS MACRO is called here. The data are saved into the file
'c:\SASWinBUGS\SAS\test\struc.txt'
*/
%_sexport(data=Sim_LinGM, file='c:\SASWinBUGS\GrowthData.txt',
var=y1-y5);
/*Create the initial values*/
DATA _NULL_;
FILE "c:\SASWinBUGS\InitValues.txt";
PUT "list(Mu=c(0,0), Inv_cov= structure(.Data = c(1,0,0,1),.Dim=c(2,2)), Inv_sig2_e=1)";
RUN;
/*Create a script file that includes the command to be used in WinBUGS*/
DATA _NULL_;
FILENAME scrip "C:\programs\WinBUGS14\bugscript.txt";
FILE scrip;
PUT // @@
#1 "display('log')"
/* Check model
1. using "/" instead of "\"
*/
#2 "check('C:/SASWinBUGS/GrowthModel.txt')"
/*Read in the data */
#3 "data('C:/SASWinBUGS/GrowthData.txt')"
/* Compile the model */
#4 "compile(1)"
/* Give the initial values */
#5 "inits(1, 'C:/SASWinBUGS/InitValues.txt')"
/* use "gen.inits()" to generate initial values for these not specified in inital value file*/
#6 "gen.inits()"
/*Discard the first 5000*/
#7 "update(2000)"
/* set variable names to be monitored, one command line for each variable*/
#8 "set(MuL)"
#9 "set(MuS)"
#10 "set(Sig2_L)"
#11 "set(Sig2_S)"
#12 "set(rho)"
#13 "set(Sig2_e)"
/*Set DIC*/
#14 "dic.set()"
/*Update to get the estimations*/
#15 "thin.updater(10)"
#16 "update(2000)"
/* call for stats to be output, one command for each variable */
#17 "dic.stats()"
/* coda file to save the coda. The data will be read into SAS later.
The index data are saved in file outputIndex.txt
The sample data are saved in file output1.txt for chain 1*/
#18 "coda(*,'C:/SASWinBUGS/output')"
/* BUGS log will be saved as bugslog.txt */
#19 "save('C:/SASWinBUGS/bugslog.txt')"
/*#20 "quit()"*/
;
RUN;
/*Creat a bat file
Can run winbugs directly. This way can save a file for further use.*/
DATA _NULL_;
FILE "c:\SASWinBUGS\run.bat";
PUT '"C:\programs\WinBUGS14\WinBUGS14.exe" /PAR bugscript.txt';
PUT 'exit';
RUN;
/*Run WinBUGS*/
DATA _NULL_;
X "c:\SASWinBUGS\run.bat";
RUN;
QUIT;
/*Read in the log file to view the DIC*/
DATA log;
INFILE "c:\SASWinBUGS\bugslog.log" TRUNCOVER ;
INPUT log \$80.;
log=translate(log," ","09"x);
RUN;
PROC PRINT DATA=log;
RUN;
/*Read in the generate samples and do the analysis
The data are saved in work.post1
*/
TITLE 'Run WinBUGS from SAS: A Growth Model Example by Zhang & McArdle';
TITLE2 'Results from the simulation study';
%coda2sas(out=post1, infile='c:\SASWinBUGS\outputIndex.txt',
chain='c:\SASWinBUGS\output1.txt', stats=1);
QUIT;
